OPEN 3 ACCESS Freely available online 



•0-PLOS I o-^E 



Quantification of Blood Flow and Topology in ^ 
Developing Vascular Networks cros^rk 

Astrid Kloosterman\ Beerend Hierck^, Jerry Westerweel\ Christian Poelma^* 

1 Laboratory for Aero & Hydrodynamics, Delft University of Technology, Delft, The Netherlands, 2 Department of Anatomy & Embryology, Leiden University Medical 
Center, Leiden, The Netherlands 

Abstract 

Since fluid dynamics plays a critical role in vascular remodeling, quantification of the hemodynamics is crucial to gain more 
insight into this complex process. Better understanding of vascular development can improve prediction of the process, and 
may eventually even be used to influence the vascular structure. In this study, a methodology to quantify hemodynamics 
and network structure of developing vascular networks is described. The hemodynamic parameters and topology are 
derived from detailed local blood flow velocities, obtained by in vivo micro-PIV measurements. The use of such detailed 
flow measurements is shown to be essential, as blood vessels with a similar diameter can have a large variation in flow rate. 
Measurements are performed in the yolk sacs of seven chicken embryos at two developmental stages between HH 13+ and 
17+. A large range of flow velocities (1 jxm/s to 1 mm/s) is measured in blood vessels with diameters in the range of 25- 
500 i^m. The quality of the data sets is investigated by verifying the flow balances in the branching points. This shows that 
the quality of the data sets of the seven embryos is comparable for all stages observed, and the data is suitable for further 
analysis with known accuracy. When comparing two subsequently characterized networks of the same embryo, vascular 
remodeling is observed in all seven networks. However, the character of remodeling in the seven embryos differs and can 
be non-intuitive, which confirms the necessity of quantification. To illustrate the potential of the data, we present a 
preliminary quantitative study of key network topology parameters and we compare these with theoretical design rules. 
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introduction 

Fluid dynamics plays a critical role in the development of 
vascular networks [1-3]. Better understanding of vascular 
development can improve prediction of this process, and may 
eventually even be used to influence the final vascular structure. 
This can for example be useful in the fields of wound healing, drug 
delivery, or tumor control. To improve insight into the relation 
between hemodynamics and vascular development, quantification 
of the hemodynamics in the developing vascular structure is 
essential. Various models of the blood vessel structure and the 
influence of hemodynamics have been described [4-6] . Validation 
of such models is crucial for further application, and this requires 
experimental data. Some quantitative information on hemody- 
namics is already available for various types of vascular networks 
(e.g. [7-1 1]), as well as methods addressing difficulties resulting 
from incomplete experimental information in networks [12]. To 
study the hemodynamic and structural changes, this data should 
preferably be available for a large network, at more than one 
moment during development. Unfortunately, such data is 
currently still lacking. This study aims to address this by presenting 
a comprehensive data set, as well as the techniques used to obtain 
it. 

The extra-embryonic vasculature of the chicken embryo (yolk 
sac vasculature or vitelline network) is a commonly-used model to 
study the development of the human cardiovascular system. 



During development, the network can be accessed multiple times 
provided that the embryo remains in the egg. This minimizes 
effects on the development of the embryo and vasculature. The 
long survival times in the in om studies reported by Korn [13] and 
the low survival rates in case of ex ovo development support this 
approach [14]. In addition, the ex om method may lead to shape 
changes of the network and thus altered mechanical stress 
distributions. As manipulation of embryo and vasculature is not 
required for the in ovo method, these detrimental effects can be 
avoided. 

The vitelline network floats on top of the yolk sac. It starts as a 
collection of blood islands [1.5]. These consist of hematopoietic 
cells (primitive blood cells), surrounded by endothelial cells, which 
will later form the vessel walls. These blood islands align and wiU 
form a primary vascular plexus, also called the capillary plexus, 
consisting of small blood vessels. As soon as the heart starts 
beating, around HH stage 10 [16], blood starts flowing through 
the early embryonic vasculature and remodeling starts [15]. 
Remodeling of the capillary plexus starts just outside the body. 
From this point, remodeling extends further outwards, and the 
capillary plexus develops into a hierarchical tree. Eventually this 
will develop into distinct arterial and venous circulations [2] . 

Hemodynamic characterization at more than one instance 
requires non-invasive measurements of the local blood flow 
velocities in a given region of interest. It has been demonstrated 
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that this can be achieved by in vivo microscopic Particle Image 
Velocimetry (micro-PIV). This technique has previously been 
applied to measure flow in the embryonic chicken heart [ 1 7 , 1 8] , as 
well as in the vitelline network [19,20]. Micro-PIV is capable of 
obtaining a velocity field in two directions (perpendicular to the 
viewing direction) on a two-dimensional grid in a planar domain 
[21]. This is sufiicient to describe the hemodynamics in the 
network: during the early developmental stages investigated in this 
study the vitelline network is predominantly confined to a thin 
planar domain. Furthermore, as the flow is laminar, the two in- 
plane velocity components wiU be sufiicient to describe the flow in 
these blood vessels. 

PrV is based on the displacements of particles in the flow, either 
artificial or naturally present. First, a series of image pairs is 
recorded to capture these particles. For the analysis, the 
measurement area is divided into smaller subdomains, called 
interrogation areas. These define the spatial resolution of the 
result. Cross-correlation is applied on the interrogation windows to 
obtain a local displacement, and eventually a local velocity. This 
statistical approach, determining the most likely displacement for a 
group of tracer particles, makes PIV a more robust method 
compared to other measurement techniques based on tracking of 
individual particles. A more extensive description of PIV can be 
found in the literature [21]. 

PIV not only provides local velocity data, but can also be used 
to obtain the network topology (see Section Methods). In principle, 
this topology could also be obtained from the brightfield images, 
which arc recorded as part of the PIV measurements. However, 
the imaging conditions sometimes do not allow a robust detection 
of smaller blood vessels, see Figure 1 . The wall location of the large 
blood vessel can be obtained by thresholding the intensity of the 
bright-field image (Figure la). However, the smaller blood vessel 
(indicated by *) is hardly distinguishable in the same figure. In 
contrast, the velocity magnitude field, as obtained by PIV, 
provides a much more distinct structure. This is illustrated in 
Figure Ic, which shows the local velocity (black dots) along the 
cross-section indicated in the figure, along with the fitted parabolic 
velocity profile (green dotted line). The corresponding pixel 
intensity values of the CCD (blue) are also shown. As can be 
seen in the figure, the vessel walls can easily be detected using the 
parabolic fitted velocity profile, but detecting the walls based on 
the pixel intensity is not straightforward. 

In this study, we make use of the capabilities of in vivo micro- 
PIV to quantify hemodynamics and determine the network 
structure in a vitelline network. Moreover, this is performed 
multiple times at developmental stages between HH 1 3+ and 1 7+ 
(2-3 day old). At HH stage 13+ remodeling has started, and 
significant changes in the network structure take place in 
progression through the subsequent stages. During the later stages 
(after HH 1 7+), the two distinct arterial and venous networks start 
to extend into three-dimensions. These stages are therefore 
excluded in this study. We focus on stages HH 13-17 and report 
data obtained from seven chicken embryos at two developmental 
stages. We evaluate the accuracy of the method and present a 
preliminary analysis of the data set to show its potential to study 
both local and global changes in the vasculature. Examples of the 
former are the fate of individual vessels, while the latter refers to a 
more statistical approach of the entire vasculature. 

Methods 

Chicken Embryos 

Fertilized White Leghorn eggs were incubated for 60 to 77 
hours at a constant temperature of 38°C and a constant humidity. 



During this incubation time, they remained in the same position, 
lying horizontally. As a consequence, the chicken embryo and the 
surrounding vitefline network were floating on top of the egg yolk. 
To perform in vivo micro-PIV measurements, optical access to the 
embryo and the vitelline network was required. For this purpose, a 
viewing window was carefully created in the egg shell [13]. The 
size of this window should approximately match the size of the 
network to enable inspection for irregularities of the entire 
network, and thus measured a few square centimeters. To prevent 
dehydration of the embryo and the vitelline network, Locke saline 
solution (0.94yo NaCl, 0.004,l"/i, KCl, 0.0()4"X, CaClj (w/v)) was 
added. To enhance the image quality, the window was then 
covered with a microscope cover glass before placing it under the 
microscope. The embryo and vitelline network were first checked 
for irregularities, and embryos with abnormal development (for 
example embryos lying on their right side or embryos with visible 
network damage) were excluded from the study. Preceding the 
velocity measurements, heart rates were measured to ensure a 
normal heart rate, and the developmental stages of the embryos 
[16] were determined (see Table 1). The stability of the heart rate 
throughout the measurement was evaluated a posteriori by 
analyzing the instantaneous velocity at the centerline of a large 
artery and counting local maxima. The heart rates during the 
measurement series are also shown in Table 1. Measurements 
using micro-PIV were performed based on bright-field imaging, as 
described in the next section. An overview image of a chicken 
embryo and surrounding vitelline network with a typical 
measurement region is shown in Figure 2. A large group of 95 
embryos were incubated for this study, but useful results were 
obtained for only 7 embryos. This small number is mainly due to 
the short time where the embryos have the desired developmental 
stages. Furthermore, the experiments required a correct position of 
the embryo and vitelline network within the shell window 
throughout the measurements. 

Following the first micro-PIV measurement s(;ri("s, the egg was 
placed back into the incubator for further development, after 
having covered the window in the egg with clear plastic tape to 
protect the embryo from dehydration. Before this, a few miUilitres 
of albumen was carefully removed through the small hole in the 
base of the egg to make sure the embryo and vitelline network 
would not stick to the plastic tape during their stay in the 
incubator. When necessary, additional Locke saline solution was 
dripped onto the egg during the remaining incubation time. For 
the second measurement series, 3 to 9 hours later (see Table 1), the 
plastic tape covering the window was removed and the egg was 
again filled up with Locke saline solution and covered with a 
microscope cover glass. 

Measurement Setup 

For the micro-PIV measurements, the egg was placed under a 
microscope in a bath containing water with a constant temper- 
ature of 38°C to ensure constant temperature and heart rate 
during the measurements (see Figure 3). An upright epifluorescent 
combi-microscope with zoom functionality (Leica MZ 16 FA) was 
used with a Leica FluoCombi III objective, providing magnifica- 
tions in the range of 2.8-46 X. A magnification of 5x with 
corresponding numerical aperture 0.14 was used. This magnifi- 
cation was chosen such that the resulting measurement resolution 
is sulficiently high while the measurement time remains acceptably 
small (up to 30 minutes for one measurement series, see the next 
section) compared to the time in which development takes place. A 
PCO Sensicam QE camera (1,376x1,040 pixels) was used for 
image recording, with 2x2 binning to improve the signal to noise 
ratio. This also enabled a two times higher acquisition rate, which 
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Figure 1. Detecting the vessel wall from both the brightfield images and the fitted velocity profile. The measured streamwise velocity 
(black dots) with fitted parabolic velocity profile (green dotted line) and the corresponding pixel intensity values of the CCD (blue) are shown for a 
typical cross section of a vessel. The pixel intensity values are averaged over 50 recorded images to improve the signal to noise ratio. The brightfield 
image and corresponding time-averaged velocity magnitude are shown on the left, with the locations of the cross sections in blue and black, 
respectively. 

doi:1 0.1 371 /journal.pone.0096856.g001 



shortened the time needed for the acquisition of the images and 
ensured tliat the pulsatile flow was adequately sampled. The 
resulting pixel size was 12.9 X 12.9 nm^, and one individual 
measurement area thus covered 1.8x1.4 mm^. The dimensions of 
the combined measurement regions are given in Table 1. 

Red blood cells (RBCs) were used as tracer, since adding 
artilicially tracer particles was less suitable for repeated measure- 
ments spanning several hours [22]. The larger size of RBCs 
(compared to typical artificial tracers) leads to an increased depth 
of correlation (DOC), which defines the measurement plane 
thickness [23]. A relatively large DOC compared to the blood 
vessel diameter can cause a severe underestimation of the 
measured centerline velocity (see [24]), since the measured velocity 
is an average of all flow velocities along the viewing direction. With 
the present magnification of 5x and the largest blood vessels 
diameter measuring 200 ^im, the measurements were in the depth- 
saturated regime [20]. In this case, the actual centerline velocity 
could easily be retrieved by applying a correction factor. 
Moreover, the measured velocities in the depth-saturated regime 
are less sensitive to small misalignments of the measurement plane 
with respect to the centerlines of the vessels compared to 
measurements using a relatively small DOC. A pulsed high- 
powered light emitting diode (LED; ILA GmbH) was used as light 
source to image the RBCs [20,25]. They appeared as black dots in 
the recorded images (see Figure 4). 

Data Acquisition and Processing 

Double-frame recording was used for the in vivo micro-PIV 
measurements using DaVis software (LaVision; Gottingen Ger- 
many), see Figure 4. This entails that subsequent image pairs were 
recorded, at a frequency of 1 / H^Iacq- Each image pair consisted of 
two recordings on separate frames, with a time difference of A/. 



The recording frequency, or repetition rate, 1 / A?^ce was set to its 
maximum, 9.7 Hz. This was independent of the heart rate of the 
chicken embryo, and hence image pairs were recorded on multiple 
and arbitrary moments during several cardiac cycles. The time 
difference between two frames of one image pair, Az, was chosen 
such that the maximum displacement occurring in all areas of one 
measurement region would not exceed 10-15 pixels during the 
cardiac cycle. As the peak velocities varied throughout the network 
and also varied with age. A? was determined for every 
measurement series separately. The maximum displacement 
restriction was achieved with A;= 10 ms, and sometimes 
A? = 5 ms or 20 ms. In some measurement regions the dynamic 
range of blood flow exceeded the measurable range, for instance if 
both very large and small blood vessels were present in the same 
field-of-view. In these cases, the velocities were also determined 
using first-frame correlation (see Figure 4). In this approach 
velocities were not based on correlating the two frames of one pair 
(ii and 12), but by correlating the first frames of two subsequent pairs 
(ii and (i+l)i). As the temporal delay between the two first frames 
is much larger than the A/ between the frames of an image pair, 
this provides a much larger dynamic range. By combining both 
correlation approaches, the dynamic range was enlarged from 
0.01-1.2 mm/s to 0.001-1.2 mm/s [22]. 

For each measurement area, 500 image pairs were recorded, 
which took less than one minute. This number was sufficient to 
apply correlation averaging [26] in order to obtain ten phase- 
locked velocity field during the cardiac cycle [22]. Correlation 
averaging based on 50 image pairs significantly improved the 
accuracy of each phase-locked velocity field, while the total 
recording time remained acceptable. After recording the 500 
image pairs, the embryo was translated using a translation stage to 
obtain adjacent regions; typically 9 areas were documented per 
embryo and stage. 
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Figure 2. The embryo and vitelline network during the third 
day after fertilization are clearly visible after windowing the 
egg. A typical measurement region is indicated by the dashed white 
box. 

doi:10.1371/journal.pone.0096856.g002 

The PIV data of multiple adjacent measurement areas, 
recorded in a single measurement series, were processed and 
combined to obtain one measurement result with the time- 





chicken embryo 
surrounded by 
vitelline network, 
on top of egg yolk 



Figure 3. For the micro-PIV measurements in the vitelline 
network, the embryo is placed in a bath of constant 
temperature of 38°C under the microscope. The vitelline network 
is clearly visible through the microscope since the network and embryo 
are floating on top of the egg yolk. A high-powered LED illuminates the 
tissue directly, not through the microscope. 
doi:10.1371/journal.pone.0096856.g003 



PLOS ONE I www.plosone.org 4 May 2014 | Volume 9 | Issue 5 | e96856 



Quantification of Blood Flow and Topology in Vascular Networks 



image pair i 



image pair 1+1 




■-ACQ 



shift correction 
phase assigning 

PIV cross-correlation 




r 



0.15 



0.1 



10.05 
0 

mm/s 



Figure 4. Subsequent image pairs are recorded to obtain a velocity field. After applying several processing steps, a two- 
dimensional velocity vector field yielding two velocity components is obtained. One measurement area is shown (embryo 3, second 
measurement series; T2), displaying the time-averaged velocity magnitude, K = |V|. A selection is enlarged to show the corresponding time-averaged 
velocity vector field, illustrating the spatial resolution of the measurement results. 
doi:1 0.1 371 /journal.pone.0096856.g004 



averaged velocity field at one developmental stage. The image 
data were processed using an in-house multi-pass PIV code [19], 
with a final interrogation area size of 8x8 pixels using 50% 
overlap. The fiaU details of this processing are described in detail 
elsewhere [22]. It entailed image drift correction, 'stitching' of the 
pulsatile velocity fields and temporal averaging. The latter was 
done here to reduce the amount of data, but transient results are 
readily available. Averaging had to be done with care at the 
appropriate step during the processing [18]: for strongly pulsatile 
flows the image data had to be sorted first based on the phase 
within the cardiac cycle; after processing, the resulting vector fields 
could be averaged. For weakly pulsatile flows the raw images could 
be processed without sorting. The velocity fields were then masked 
to exclude regions with outliers (generally occurring outside the 
blood vessels) by a median filter. The end result of the PIV analysis 
and post-processing was, for each embryo and stage, a time- 
averaged velocity field of typically 3x5 mm^ with a relatively high 
spatial resolution of 10x10 fim^. 



Hemodynamic Parameter Extraction 

To facilitate a quantitative analysis, each measured two- 
dimensional velocity field was modeled as a coUection of 
connected vessel segments. Vessel segments were each character- 
ized by a single diameter, velocity and length. Diameter and 
velocity were obtained by averaging results of multiple adjacent 
extracted blood vessel velocity profiles. Based on these parameters, 
other relevant hemodynamic parameters can be calculated, such 
as e.g. the flow rate and the wall shear stress. Typical values for 
embryonic avian blood properties, required for the determination 
of the wall shear stress, are reported in the literature [27]. 

To define the network structure, blood vessel branches were first 
automatically identified based on the skeleton [28,29]. Vessel 
segments between two branch points were subsequently labeled 
using a Matlab routine [22,30]. As an illustration, a measurement 
region with corresponding time-averaged velocity magnitudes is 
shown in Figure 5. The non-linear color-coding was adapted such 
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that all velocity scales are visualized. The skeleton is displayed as a 
black line, branch points are marked by white circles, and end 
points are marked by gray circles. Accordingly, the skeleton's part 
between two nearest branch points (including end points) 
approximately coincides with the centerline of one vessel segment. 
This way, the collection of vessel segments and their connectivity 
was defined. 

Next, a characteristic velocity and diameter for each vessel 
segment were extracted. The procedure is explained using one 
vessel segment, shown in the top part of Figure 6. An extensive 
description of the procedure is given by Kloosterman [22]. To 
obtain a single characteristic vessel diameter and velocity, the flow 
rates were assumed to remain constant throughout each vessel 
segment. The entrance length of the flow can be considered to be 
small compared to the length of a vessel segment due to the low 
Reynolds number (Re<<l), and the flow could therefore be 
considered fully-developed throughout each vessel segment. A 
polynomial fit was used to obtain vessel diameter and centerline 
velocity. At least three points are needed for a polynomial fit. This 
limited our approach to blood vessels larger than three times the 
measurement resolution, i.e. approximately 30 fim. 

Poelma et al. (2010) [18] showed that the approach to 
determine the wall location based on the velocity profile could 
be accurate within a few micrometers in a comparable in vitro 
experiment. Avian blood has a relatively constant viscosity, in 
contrast to the well-known strongly shear-thinning behavior of 
human blood [27]. This implies that the velocity profile wfll be 
nearly parabolic, which has indeed been confirmed in previous 
experiments [20] . At multiple cross sections of a vessel segment 
(white and red lines in the top graph), parabolas (second-order 
polynomial) were fitted to the local velocity profiles, which gave 
local centerhne velocities (Fcl) and diameters d. Neighbourhoods 
of branch points were omitted in this procedure since defining the 
streamwise direction here is not straightforward. The trend 
showed that the centerline velocity (Fcl, on the bottom of the 
graph) was inversely proportional to the cross-sectional surface 
area (A~d^, middle part of the graph). An average flow rate Q, 
was determined based on the valid flow rates (flow rates between 
0.5 and 1 .5 x the median flow rate). These data points are located 
in the figure in the green valid region; excluded data points are 
encircled red. The standard deviation of the valid flow rates was 




mm/s 



Figure 5. Time-averaged velocity magnitudes of a measure- 
ment region are sKiown, together with with the corresponding 
skeleton (blacl< line), branch points (white circles), and end 
points (gray circles). The flow enters the measurement region from 
the right. Note that a non-linear color scale was used for the velocity 
magnitude. 

doi:1 0.1 371/journal.pone.0096856.g005 



generally an order of magnitude smaller than the average flow 
rate. Now, the characteristic centerline velocity was determined by 
averaging the velocities from the valid data points. The 
characteristic diameter was subsequently derived from the average 
flow rate and characteristic centerline velocity. Additionally, the 
result was visually inspected to detect any remaining anomalies 
that might still be present. 

Data Quality Assessment 

It was not possible to directiy validate the measurement 
technique in vivo. In an in vitro validation experiment, the 
technique has previously been shown to be capable of measuring 
flow rates with an error less than 2% error [18]. For the current 
in vivo application, the only way to assess the accuracy of the 
velocity measurements and the subsequent segment identification 
was to perform internal checks. An example is the aforementioned 
variation of the flow rate in a single vessel, which can only be due 
to measurement errors. To compare the performance of the 
parameter extraction method, the cross-sectional velocities and 
corresponding fitted velocity profiles for three vessel segments were 
investigated in detail. These segments were selected to have a 
similar diameter, but with significantly different velocities. 

The quality of the total data set was also investigated by 
verifying the flow balances in all branch points [31]. The flow 
balance is defined by Qin — Qom, and the relative flow balance by 
(2/1! ~ Qoiii)l Qin, where Qi„ and Qom is the total flow entering and 
leaving the branch point, respectively. An ideal flow balance 
without measurement errors is obviously established when 
Qin = Qout- For every measurement region the distribution of flow 
balances was studied. A relative small number of large outiiers, 
where |(g,„ - Qout)l Qin\> 1, was excluded. 

Our method is based on representing the flow based on a 
polynomial fit. Measurement errors are more likely to occur near 
the wall, where the velocity gradients are at a maximum. 
Furthermore, fits in larger vessels are based on more velocity 
information. The uncertainty in the fit - and thus flow rate - can 
therefore be expected to decrease for increasing vessel diameters. 
To check this, we evaluated the flow balances for both the fuU 
collection of blood vessels, as well as a subset containing only the 
larger vessels. These were identified as vessels larger than 52 /im (5 
times the spatial resolution of the velocity measurements). This 
subset was roughly three-quarters of the whole set of vessels. 

Structural Parameters 

In this study, six additional parameters were calculated to study 
structural changes between two consecutive networks. The first 
two parameters were related to the total network structure: the 
vessel occupancy and the vessel density. These were defined as the 
part of a measurement region occupied by vessel segments, and 
the average number of vessel segments present in one square 
millimeter. The tortuosity T, defined as the ratio of the actual path 
length L and the Euchdean distance between the two adjacent 
nodes Lo, as well as the angular vessel flow 9, defmed as the angle 
between the main flow direction in the total measurement region 
and the local vessel flow direction (see Figure 7 for illustrations), 
were determined for all vessel segments. 

Also the bifurcations were inspected. Parameters derived from 
design principles for optimal networks were validated by the 
experimental values from this study. Every branch point connects 
one parent vessel with two child vessels having the same flow 
direction: towards (converging, in this study defined as venous) or 
away from (diverging, in this study defmed as arterial) the branch 
point. Note that the arterial and venous bifurcations were defined 
in this study by their flow directions (diverging or converging. 
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Figure 6. Multiple parabolic fits at several locations on the vessel centerline were performed to obtain a single characteristic 
velocity and diameter for each vessel segment. The time-averaged flow rate is assumed constant throughout the vessel segment. The valid 
region is bounded by 0.5 and I.Sxthe median flow rate, and the red-encircled data points were not incorporated, due to a strongly deviating flow 
rate. Note that the fitted diameters and flow rates for the two data points on the far right are too large to be visible in the graph. 
doi:1 0.1 371/journal.pone.0096856.g006 



respectively), and not by specific gene expression. First, the 
branching angle a was determined. This was defined by the angle 
between the two lines connecting the branch point with 
respectively the center of one daughter vessel halfway the two 
adjacent branch points, and the center of the other daughter vessel 
halfway the two adjacent branch points (see Figure 7 for schematic 
diagrams). 



Much research has been done on the validation of Murray's law 
([32-34], among others), also in the vitelline network of the 
chicken embryo [35,36]. The design principle of a minimal cost 
function leading to this law was first introduced by Murray in 1926 
[37,38], and holds the relation between the parent vessel diameter 
-Do and child vessel diameters Di and D2: -Dq = Z>j+Z)2- To 
vaUdate this law for our data sets, the parameter y = Z)q/(Z)j -I--D2) 
was computed, and should ideally be equal to one. Since this 




Figure 7. Schematic diagrams of studied structural parameters are shown for clarity. These parameters are: the tortuosity t = L/Lq, the 
angle between the vessel segment and the main flow direction <Jo, the arterial branching angle a„rt_„,,g, the venous branching angle o!„,„_„,j;, and 
Murray's law ratio for both arterial and venous branches }'„,.,j„,g and }'„,„_„,,g. 
doi:1 0.1 371/journal.pone.0096856.g007 
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parameter is veiy sensitive to outliers, tlie impact of outliers was 
reduced by removing the largest outlier from the data set ()' 
exceeding the value of three). The value of three was high enough 
to exclude the large outliers affecting the statistics; lowering this 
limit had a minimal effect since only a few outliers were removed. 

Results and Discussion 

The vitelline networks of the seven embryos at two consecutive 
developmental stages (Tl and T2) included in the study were 
modeled as a collection of connected vessel segments, each with 
one characteristic diameter, velocity, and length. The complete 
collection of modeled networks of aU seven embryos can be found 
in Appendix B in Ref [22]. The data (File SI) and visualizations 
(Figure SI) of all networks are also available as electronic 
supplementary material of this manuscript. 

Data Quality 

For the three vessel segments A, B, and C in Figure 8, the cross- 
sectional measured velocities and corresponding fitted velocity 
profiles are shown in Figure 9. The diameters are in the same 
range: between 106 and 131 /urn. The velocity differences, 
however, are considerably larger: the velocity in vessel segment 
A is almost ten times the velocity in vessel segment C, but a clear 
parabolic shape of the cross-sectional velocity is observed for all 
three vessel segments. The comparable performance of the 
method in a large velocity range clearly indicates that the variation 
in the extracted velocity in one measurement region is caused by 
biological variance in dimensions and flow, and is not a 
consequence of inaccurate measurements or fitting procedures. 
This also implies that estimating local blood velocities from vessel 
diameter (for example by Pries et al. (1990) [39]) may lead to 
errors; as can be seen in Figure 9, a vessel with a diameter of 
1 00 fjm can have a centerline velocity that varies over an order of 
magnitude. 

The distribution of the relative flow balance of each measure- 
ment region could be approximated by a normal distribution, 
when excluding a relative small number of large oudiers where 
KQin — Qout)/ Qm\> ^- The bar chart in Figure 10 shows the 
probability density function of the relative flow balance for embryo 
3 at T2; the thick black curve represents the fitted normal 
distribution. Additionally, the corresponding mean ii and bound- 
aries of the respectively 68% and 95% tolerance intervals (/i + c 
and jJ. + 2a) are indicated. According to these parameters, 68% of 
the branch points have a relative flow balance between —30% and 
40% for this measurement region. A summary of the mean, 
standard deviations, and percentage of outliers can be found in 
Table 2 (listed under 'all diameters'). The mean values are close to 
the theoretical value of zero for all embryos, so there appears to be 
no systematic flow measurement error. Comparing the statistics 
for the whole set and the 'large diameter' subset, there seems to be 
no difference. This suggests that the relatively larger errors in 
smaller vessels (as anticipated in the Section Data Quality Assessment) 
appear to be absent. Note that the standard deviations of the 
relative flow balance a reported are significantly larger than the 
standard deviations of the flows that are used to calculate the 
balance. Error progression can easily double the error in the 
reported values for a. This means that a conservative estimate of 
the standard deviation of the flow Q,is approximately 20%; this is 
somewhat higher than the variation along the centerline reported 
earlier (see e.g. Figure 6). 

The two cases that deviate most from the theoretical situation 
and also have the largest relative numbers of outliers are the 
measurements of embryo 1 at Tl and embryo 7 at T2. For 




Figure 8. Networks are modeled as a collection of connected 
vessel segments, here shown with corresponding diameter and 
color-coded time-averaged mean velocity for two consecutive 
measurements T1(a) and T2(b). Besides a general velocity increase, 
the part of the network in the upper right corner seems to have 
become responsible for a larger part of the total blood flow in this 
region (see also text). The flow enters the measurement region from the 
right. Note the two different non-linear color scales for (a) and (b). 
doi:10.1371/journal.pone.0096856.g008 

embryo 1, this can be explained by the fact that at the time of 
measurement, this embryos had only reached HH stage 13-H. At 
this stage, the vitelline network stiU resembles a sponge-like 
structure [2,40] and it is difficult to define individual blood vessels 
and thus network topology. For the case T2 of embryo 7, the 
embryo had reached HH stage 17-H. At this stage, the network 
starts to extend into three-dimensions instead of being confined to 
a plane. This is due to the growth on the curved yolk sac. The 
measured velocities are hence a projection of the actual velocities 
onto the measurement plane, and are likely to be slightiy 
underestimated. This also results in a less accurate overall flow 
balance. 

Example of Network Remodeling 

Figure 8 shows the two representations for embryo 3. This 
embryo developed from HH stage 14 to 16 during this study, 
which took approximately 4 hours. Despite this apparently short 
time, significant differences in flow and structural changes 
occurred. The maximum time-averaged velocity in this measure- 
ment region more than doubled (with an even greater increase in 
volumetric flow rate). The preferred direction also seems to have 
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Figure 9. The variety in time-averaged mean velocities for all 
vessel segment diameters are shown for two consecutive 
measurements T1 and T2. For three vessel segments (I, II, and III, 
indicated in Figure 8), the changes in diameter and time-averaged 
velocity from T1 tot T2 are indicated by the dashed arrows. The cross- 
sectional velocity profiles shown for three data points (A, B, and C, 
indicated in Figure 8) show consistency despite the different time- 
averaged velocities. 
doi:1 0.1 371 /journal.pone.0096856.g009 



changed: the part of the network in the upper right corner seems 
to have become responsible for a larger portion of the total blood 
flow. For this embryo, it means that the preferred direction of the 
blood flow shifted towards the tail of the embryo. The changing 
hemodynamics are also clearly visible in the upper part of Figure 9, 
where time-averaged mean (bulk) velocities are shown as function 
of the respective vessel segment diameter for the two stages Tl and 
T2. For three vessel segments (I, II, and III, indicated in the top 
part of Figure 8), changes in diameter and time-averaged velocity 
from Tl tot T2 are indicated by the dashed arrows. The 
development of these three vessel segments show completely 
different behavior. Where vessel segment I splits into two vessel 
segments due to a new branch midway along the segment, the 
time-averaged velocity and flow in vessel segment II more than 
doubled with minimal change in diameter. On the contrary, the 
time-averaged velocity and flow rate in vessel segment III 
decreased. 

Qualitative and Quantitative Observations 

For every embryo in this study, changes in structure and 
hemodynamics took place between the first and second measure- 



o 

E 

a; 



5 
o 



n 



E 



? E 



^ I 
at .2 
u -a 

<u S 

E 



T3 

V 5 



c c 



E »! 



to p 
fit ^ 



E 



° E 



I Ch OJ 



Q> fU OJ 

-2 ^ , O 

re qj c: O 

o fD in g 

^ o 

OJ OJ -M o 



o Oi t 



Ol .9 J! 



PLOS ONE I www.plosone.org 



9 



May 2014 | Volume 9 | Issue 5 | e96856 



Quantification of Blood Flow and Topology in Vascular Networks 




Figure 10. The bar chart shows the probability density function 
of the relative flow balance (Q/,, — Q„iii)IQm for embryo 3, at T2. 

This can be represented by a normal distribution which should ideally 
have a mean /i close to zero and a small standard deviation a. 
doi:1 0.1 371 /journal.pone.0096856.g01 0 

ment series. Direct comparison between the embryos was difficult, 
since the measurement regions of the seven embryos were not 
always located in the same part of the vitelline network. However, 
one could conclude that remodeling took place in all develop- 
mental stages in this study, and a rich array of observed 
phenomena was observed. The qualitative changes were com- 
pared with the quantitative results, and are described in the 
following paragraphs. 

The distributions of the characteristic diameters D, time- 
averaged velocities V, and lengths L of the vessel segments were 
represented by the 25%-50%-75% percentiles for all measure- 
ment series (see Table 3). During development, minimal changes 
in the percentiles of the diameters were observed (up to 23% 
dilference in the percentiles). However, the overall velocity 
increased by a factor of three in some cases. This large velocity 
increase and a minimal heart rate increase were measured for 
embryo 7. Despite the fact that the mean velocity can be 
influenced by the heart rate, the changes in velocity due to 
development are found to be more substantial. 



During development, a network can remodel into a denser 
network. This was the case for embryo 5: the network at T2 could 
be described as the network at Tl with additional smaller vessel 
segments filling up areas between the existing vessel segments. An 
increase in the number of vessel segments and vessel density (see 
Table 1) were evident consequences. These observations were also 
consistent with the changes in other parameter distributions: the 
measured diameter and velocity distributions showed an increased 
presence of small blood vessels with corresponding lower velocities, 
and also an overall shift to shorter vessel segments was measured. 
This overall decreased vessel length was caused by longer vessel 
lengths breaking up into multiple smaller vessels due to new 
branches, and obviously, the new small vessels themselves. A 
relatively large increase of the vessel density with respect to the 
vessel occupancy also indicated presence of smaller vessels. 
Furthermore, these new vessel segments connected the previously 
existing vessel segments, hence they will have a tendency to be 
directed away from the main flow direction. The resulting 
increased spreading in vessel flow direction can be related to the 
greater standard deviation of the vessel flow direction ag (see 
Table 4). Embryos 3, 6, and 7 had an approximately similar 
progression, with the primary difiFerence being an increase in the 
overall velocities. This resulted from greatiy increasing velocity 
magnitudes in the larger vessels. From a physicogical standpoint, a 
denser network may be necessary to expand the exchange area for 
nutrition supply, or, at a later stage, part of these smaller vessels 
may form the venous network [2] . 

In contrast to the above general trends, a network can also 
remodel into a less dense network. This is likely to happen earlier 
in development: vessels disconnect or merge to make the network 
more efficient for transport. This was the case for embryos 1 
(HH13+ HH14+) and 4 (HH14 HH16). Similarly, the 
measured hemodynamic parameters indicated this process by a 
decrease of vessel density and even vessel occupancy, and an 
increase in the overall vessel lengths. 

Development and Structural Parameters 

Table 4 lists statistics of several parameters for all measurement 
series, including vessel tortuosity (t), flow direction [6), and 
branching angle (a). While visually the tortuosity seemed to have 
changed during development, the statistics of the measured 
tortuosities did not confirm this. These structural changes could 
probably better be described as smoothing of the vessel segments 
with decreasing curvatures. Nevertheless, one notable feature of 
these statistics was the almost constant mean and standard 



Table 3. Distributions of the characteristic hemodynamic parameters for every embryo. 





embryo 


percentiles D (/im) 




percentiles K(//m/s) 


percentiles L (mm) 




label (no.) 


Tl 


T2 


Tl 


T2 


Tl 


T2 


1 


56-82-119 


52-76-119 


25-46-103 


18-36-123 


0.12-0.18-0.27 


0.11-0.20-0.33 


2 


52-81-115 


46-66-91 


14-43-63 


25-44-87 


0.14-0.23-0.33 


0.14-0.24-0.36 


3 


50-73-112 


48-64-93 


20-35-55 


21-33-62 


0.11-0.22-0.45 


0.12-0.20-0.32 


4 


51-77-97 


50-74-97 


33-58-92 


18-32-50 


0.14-0.22-0.37 


0.14-0.25-0.39 


5 


57-82-98 


45-71-105 


23-61-106 


13-33-77 


0.13-0.22-0.40 


0.10-0.16-0.27 


6 


62-79-99 


63-82-107 


32-53-83 


66-119-200 


0.09-0.16-0.25 


0.09-0.15-0.22 


7 


54-74-97 


45-63-83 


33-61-94 


103-192-321 


0.09-0.15-0.23 


0.09-0.14-0.23 



For every embryo, the distributions of the characteristic diameter D, characteristic time-averaged velocity V, and characteristic length L of the vessel segments are 
respresented by the 25%-50%-75% percentiles for both measurement series 11 and T2. 
doi:l 0.1 371/journal.pone.0096856.t003 
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^ £ 



O c 



^ _g deviation for all measurement series. Regardless of the stage, the 

J average tortuosity was close to 1.1, which indicates that an extra 

^ I 10% vessel length is added to the shortest path to obtain the actual 

g -5" vessel length. 

.| The standard deviations of the mean branching angle for each 

-o 5 set 0£ varied between 15 and 38°, with no clear distinction between 

^ 'j; arterial or venous junctions. The changes in branching angles 

g "5 during development also showed no clear trend, and no 

2 considerable differences for this parameter between the arterial 

^ S and venous branching points in the same measurement region 

II g were observed. However, the optimal branching angle, based on 

-a the theory for both optimal symmetrical and nonsymmetrical 

o — branching should be approximately 90° [41]. This value was 

o closely approximated by this data set. Besides, assuming that a 

^ _^ branching angle could have any value between 0° and 180°, with 

■q E probabihties equally distributed for all possible outcomes, the 

c £ standard deviation would be 52° (180/ \/T2, deduced for a uniform 

'% distribution). The standard deviations for the different data sets 

^ S ranged from 15-38°, with an average of 25°. This indicates that 

"E y> the branching angles were not completely random, which would 

? t have resulted in a mean and standard deviation of 90°±52°, but 

TOD 

S were aligned with the flow direction. 

The standard deviations corresponding to the mean values of y 
varied between 0.4 and 1.2 (except for the second measurement 
region for embryo 7). This was more than three times the 
difference in mean values of y for arterial and venous branches for 
each measurement region. In addition, the variations between the 
different measurement regions were found to be minimal, with the 
exception of the second measurement region for embryo 7. This 
could be caused by the decreased validity of the two-dimensional 
approximation of the planar vitelline network. The network 
extends over a larger part of the curved yolk sac due to growth. 
The determination of -D will be affected by these curvature effects. 
Since is involved in this law, small errors in diameter extraction 
have a large effect on the value for y. Nevertheless, the theoretical 
value of unity was in most cases closely approximated by the 
average value of this parameter. No considerable diflFerences for 
this parameter between the arterial and venous branching points 
in the same measurement region were observed. 



Conclusions 



c !s A methodology is described where in vivo micro-PIV measure- 

I, e ments are used to quantify hemodynamics and network structure 

^ g in relatively large vascular networks. Since the measurement 

S £ method allows repeated measurements, it can be applied in a 



network multiple times during development. To demonstrate the 



.2 possibilities of the method, it has been applied on the yolk sacs of 

^ seven chicken embryos at two developmental stages. The 

"g J developmental stages are in the range of HH stage 1 3+ to 1 7+, 

^ ■§ after remodeling has started and before arteriovenous diflferenti- 

^ 'I £ ^ ation has occurred. The measured regions of the vitelline network 

"g _g S had an average surface of 13 mm^. In total, 3901 vessel segments 

= -D g, S were characterized. The measured vessel diameters were in the 

S I! ^ g range of 25-500 ^im, and the measured time-averaged mean 

■K £ "g ? velocity in the range of 1 fim/ s - 1 mm/ s. With this measurement 

S S" S o technique, a large range of flow rates spanning four orders of 

o" = 1 magnitude were observed. A large variation of flow rates (order of 

J~ > ^ g magnitude) can be present for diameters in the same diameter 

§ :5 c ^ range. This exemplifies that local flow rate measurements (e.g. by 

I" I S g micro-PIV) can provide additional information, and can have 

g I E great benefits compared to studies solely based on imaging. 
Parameters such as the flow rate and wall shear stress can easily be 
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calculated from the measured parameters, for example to 
investigate vascular development. 

To verily the quality of the data sets, the accuracy of the 
extracted characteristic vessel diameter and time-averaged velocity 
were determined by evaluating the relative flow balances in the 
branch points. Ideally, these relative flow balances should be equal 
to zero, based on mass conservation. Generally, the relative flow 
balance distribution could be approximated by a normal 
distribution around zero, when large outliers were excluded. This 
analysis showed comparable results for all measurement regions, 
and an estimated standard deviation of the flow rate of 
approximately 20%. Furthermore, the two networks that attract 
attention are the network during early remodeling (HH13+), and 
the oldest network in this study (HH17+). The network in HH 
stage 1 3+ is still a sponge-like plexus instead of a tree-structured 
network. This complicates defining individual blood vessels, which 
is likely to influence the relative flow balance. The network in HH 
stage 17+ starts to extend into th(" third dinK'nsion, which is also 
likely to influence the relative flow balance. Nevertheless, the 
quality of the data sets from the seven embryos is comparable, and 
the data is suitable for further analysis with known accuracy. 

We performed a preliminary analysis to highlight the potential 
of the data set for further research on development. The 
qualitative observations were compared with the changes in the 
quantitative parameters. Vascular remodeling was observed in all 
obtained data sets, but the character of remodeling differed. For 
example, while some networks became denser, other networks 
remodeled into a less dense network. Likely, the dilferent locations 
in the total \itelline network and the developmental stages 
influence the character of remodeling. Despite the dilferent 
behavior, the changes in network parameters, such as the vessel 
diameter, flow velocity, length, and flow direction could be related 
to the observed structural changes. The agreement of the data sets 
to design rules for obtaining optimal networks were also 
investigated. For this, the branching angle and the agreement 
with Murray's law were evaluated. For all individual data sets, 
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